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A variant of White's density matrix renormalisation group scheme which is designed to compute 
low-lying energies of one-dimensional quantum lattice models with a large number of degrees of 
freedom per site is described. The method is tested on two exactly solvable models — the spin-1/2 
antiferromagnetic Heisenberg chain and a dimerised XY spin chain. To illustrate the potential of 
the method, it is applied to a model of spins interacting with quantum phonons. It is shown that 
the method accurately resolves a number of energy gaps on periodic rings which are sufficiently 
large to afford an accurate investigation of critical properties via the use of finite-size scaling 

theory. 



I. INTRODUCTION 

Since its inception, White's density matrix renormali- 
sation group (DMRG) methodlil has proven to be a pow- 
erful, robust and portable numerical method for calcu- 
lating properties such as excitation energies and correla- 
tion functions of low-dimensionalr«uantum lattice mod- 
els in condensed matter physicsEI. The first applica- 
tions of the method were to the calculation of static, 
zero temperature properties^^one-dimensional models 
with short range interactionstffl'ffiffl, but the-onethod has 
been extended to include disordered systemscl, dynamical 
propertiestj, classical two-dimensional systemsEE —finite 
temperature propertiesEil, non-hermitian systems£3, and 
systems with long range interactions^. Further-mace , ni-. 
oneering extensions to two spatial dimensions^OEjO'Ej 
have rendered the method competitive with any other 
method for key-two-dimensional models such as the Hub- 
bard modelMll3. The domain of applicability of the 
method continues to. grow with recent applications to 
wetting phenomenaEy and quantum chemistryEil. 

In principle, little modification is needed in order to ap- 
ply the DMRG to models with a large number of degrees 
of freedom per lattice site, such as models with bosonic 
degrees of freedom, an example being electron-phonon 
models. As with exact diagonalisation calculations on 
small clusters, it is possible to truncate the single-site 
Hilbert space (e.g., by limiting the phonon number in 
an electron-phonon model) to the extent where the error 
in ignoring the shed states is negligible. In practice this 
will generally mean that the "added sites" in the DMRG 
calculation!!! will contain considerably more states than 
for a spin or fermion model, but calculations have been 
performed in this manner for the Bose-Hubbard modelE3, 
a model involving spins interacting with phononsEj, as 
well as pure acoustic phononscj. For some models, how- 
ever, the number of states required to accurately repre- 
sent a single site can begin to become comparable to the 



number used in representing a whole block (i.e. the "sys- 
tem" or "environment" blocka), and this straightforward 
implementation of the DMRG can become inefficient or 
unworkable. 

White and co-workers have developed two methods for 
dealing with the single-site Hilbert space. In one ap- 
proach, 2™ truncated phonon Hilbert space levels are rep- 
resented by u spinless fermions and a standard "sweep- 
ing" methodu is employed to add just one fermion site 
at a time rather than a whole phonon degree of free- 
dom. This method was used to study the polaron prob- 
lem (a single electron interacting with a lattice of vi- 
bratins-atoms) in the one- and two-dimensional Holstein 
modelsE3. A more promising approach, that of "local 
Hilbert space reduction"t3, involves finding a highly ef- 
ficient single-site basis which can be used as an alterna- 
tive to simple Hilbert space truncation when performing 
DMRG or exact diagonalisation studies. This is achieved 
by taking a small system that can readily be exactly 
diagonalised, e.g. a system with 4 sites, and then us- 
ing the projection operator for the ground state and/or 
some low-lying excitations to define a reduced density 
matrix for a single site, by tracing the projection op- 
eratef, over the degrees of freedom of all but one of the 
sitesca. The resulting basis was shown to be very efficient 
in exact diagonalisation studies of four- and six-site half- 
filled Hubbard-Holstein systems, in that the number of 
states-required to accurately represent the site was very 
smaller. 

In this paper we present a DMRG scheme for one- 
dimensional quantum lattice systems with a large number 
of degrees of freedom per site which has some-similarities 
to the local Hilbert space reduction schemeO, and also 
to Wilson's computer renormalisation group methodLj. 
We will call the method the "four-block method" , as 
it uses four blocks rather than the two blocks used in 
the standard DMRG. The four-block method has been 
used to perform a study of an electron-phonon model — 
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the spinless fermion Holstein modelEa. In this study the 
phase boundary separating the metallic and insulating 
phases of the model was determined with high accuracy, 
and very good agreement was obtained with analytical 
results which become exact in the nontrivial strong cou- 
pling limitEl The four-block method is described in Sec- 
tion H The accuracy is then tested on solvable models. 
Finally, the potential of the method to study systems 
with a large number of degrees of freedom per site is 
demonstrated by presenting some convergence results for 
the XY spin-Peierls model. 



II. THE FOUR-BLOCK METHOD 

The four-block method is illustrated schematically in 
Fig. [l[ A calculation commences with a ring of four sites. 
A basis for the ring is the product of 4 copies of a single- 
site basis i.e. 



® B, <g> B, 



B a 



(1) 



where B s = {|n) : n = 1, 2, 3, . . .} is a single-site basis. 
For example, for a pure phonon system we might choose 
|n) to be the nth phonon level. The ring basis is reduced 
down to a finite set by limiting the size of the single-site 
basis. That is, B llng is replaced by 

S ring = B s m) ® 5 s (m) <8 # s (m) <8 # s (m) (2) 

= {\ni) <g) |n 2 ) <8> |n 3 ) (g> |n 4 ) : 

ni,n 2 ,n 3 ,n 4 = 1, . . . , m} (3) 

where B^ = {|n} : n = 1, . . . , m} is a truncated single- 
site Hilbert space. The cutoff m is chosen so that, for 
the purpose of finding low-lying excitations of the ring 
Hamiltonian, the error in ignoring the shed states is neg- 
ligible. The ring Hamiltonian is diagonalised by a sparse 
matrix method to produce the ground state |^). 

The next step is to consider a "system" block A con- 
sisting of sites 1 and 2, as shown in Fig. [I]. A truncated 
basis for this block is 



A = B ring 
= {\ni) 



r>( m ) 
ring 

|n 2 ) : ni,n 2 



1. 



,,m}. 



(4) 
(5) 



A reduced density matrix pa is defined for block A by 
integrating out the "environmental" degrees of freedom, 
n 3 and n 4 , from the ground state projection operator 



(nil <g> (n' 2 \p A \ni) <g> |n 2 ) 

m 



/ j rn\n 2 ri3ri4 

71 3 ,714 — 1 



(6) 



where ^ ni „ 2 „ 3 „ 4 = {(rc.i| ® (n 2 \ <8> (n 3 | ® (n 4 |} |^>). Next, 
Pa is diagonalised yielding eigenvalues cj„ and eigenvec- 
tors |n}}. The density matrix eigenvalues w„ are real and 



positive and sum to unity, as pA is a probability matrix, 



1 > u>i > u>2 > • • • > w m 2 > 0; 

m 2 

= 1- 



(7) 
(8) 



At this stage a cutoff fa < m 2 is chosen and a new trun- 
cated basis B^ is developed for the block A from the 
m "most important" density matrix eigenstates, i.e. 



Bf)={|n»:n=l ) 



(9) 



The procedure of forming a ring, or "superblock" , is 
then repeated, this time using four two-site blocks (copies 
of the system block) instead of four single sites, as shown 
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B^ and |n)) play the roles of m, B^ and 
|n), and again the superblock Hamiltonian is diagonalised 
for the ground state |^); a new system block, consisting of 
two blocks (or 4 sites) is generated, and a reduced density 
matrix is defined. The procedure is iterated, with the 
lattice (or superblock) size doubling at each iteration, as 
is the case with Wilson's computer renormalisation group 
methodc3. The method also has some similarity with the 
local Hilbert space reduction techniques in that initially 
a local Hilbert space is defined for a two-site block by 
means of a reduced density matrix. 

At first glance the four-block method might appear dif- 
ficult to implement as the size of the superblock Hilbert 
space grows (with the number of states retained per 
block, m) as m 4 rather than to 2 for the standard DMRG 
methodu. Fortunately, there are a number of steps that 
can be taken to reduce the computer resource require- 
ments. 

1. In defining a truncated basis for the superblock, 
rather than choosing m and taking 



(m) 



n( m ) = 

super 



{K» ® M) ® K» ® |n 4 » : 
ni,n 2 ,n 3 ,n 4 = l,...,m}. 



(10) 



we choose a cutoff < e < 1 and define the trun- 
cated superblock basis according to 

S s ( 4er = {K» ® K» ® |n 8 » 8> |n 4 )) : 



>^}. (11) 



That is, the superblock basis is taken to be the 
set of all fourfold products of density matrix eigen- 
states such that the product of the corresponding 
density matrix eigenvalues is e or greater. Thus, e 
is the single truncation parameter which determines 
the accuracy of the calculation. Denoting the size 

of the superblock Hilbert space by = i?i u per , 

we note that or the accuracy of the calcula- 

tion, increases as e is decreased. Typical values of 
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e range from 10 -10 to 10 -25 , and allow the study 
of systems with very large effective values of to, i.e. 
some hundreds of states can be retained per block. 
For a given accuracy requirement, the advantage 
(in terms of CPU and memory) of using £?iupcr in- 
stead of i?iupcr is some orders of magnitude. This 

is because many of the states in -Biupcr contain two 
or more block states with low density matrix eigen- 
values and thus have very low probability. These 
states are not considered when fiiuper is used. This 
approach to Hilbert space truncation can also be 
applied in the standard DMRG algorithm. A mod- 
ification along these lines is given irO. 

2. Because the four blocks that make up the su- 
perblock are identical, use can be made of the 
translational and/or reflection symmetries in re- 
ducing the CPU time required to act the superblock 
Hamiltonian 7i S upcr on a state. This operation, 
needed for the sparse diagonalisation of 7i SU pcr, i s 
the most CPU intensive operation in the algorithm. 

3. The algorithm has a natural vectorisation. Be- 
cause any given term in 7i sup cr only connects two 
of the four blocks, inner loops can be taken over 
the states of the inactive blocks. For example, 
consider a term in W supcr , 7ii2, which connects 
blocks 1 and 2. Taking an initial superblock state 
l^i) €5 1«2) ® \na) €5 1^4), the most general final state 
(under the action of H12) is l^i) ® \n' 2 ) ® |ri3) <8> \ri4), 
i.e. the indices for blocks 3 and 4 are unchanged. 
Making outer loops over rii, n 2 , n[ and n' 2l the 
matrix element \ = ( n 'i I® ( n 2 \T~Li2 1 n i) ®| n 2 ) is 
calculated or read in from storage. An inner loop 
in which \ is repeatedly reused can be performed 
over 713 and 714 (melded into a single index). The 
superblock states can be ordered in such a way 
that all memory access (to arrays representing su- 
perblock states) in the inner loops is contiguous. 



III. ACCURACY TESTS FOR THE HEISENBERG 
CHAIN 

The four-block method has been tested on the 5=1/2 
antifcrromagnetic Heisenberg spin chain 



JV 



t~L — 2 Si.Si+i, 



(12) 



where Si is the spin-1/2 operator for site i and a periodic 
ring of N sites is .assumed. This model is exactly solvable 
by Bethe ansatzed and in particular, exact results are 
available for the ground state energy Eqs , and the,singlet 
A ss and triplet A st gaps on finite periodic ringsEx 

The states associated with these gaps are found by 
making use of the parity (spin-flip) operator 



f : 5f > Sf , S i > 5j , (13) 
and utilising a projection operator of the formO 



instead of the ground state projection operator to define 



the density matrix, where 



(0) 



4 1} 



and 



4 0) \ de- 



note the ground state and the first excited triplet and 
singlet states respectively. 

Results for Eqs , A ss and A st for various values of e and 
N = 32, 64 and 128 are given in Table |. The high accu- 
racy for N — 32 is to be expected as a substantial frac- 
tion of the complete Hilbert space is retained in this case. 
For N — 128 the singlet and triplet gaps are resolved to 
within around 0.1%. Such accuracy for periodic rings in 
critical models where the gaps vanish in the thermody- 
namic limit makes the four-block method potentially use- 
ful for finite-size scaling studies. The four-block method, 
like the standard DMRG, is variational in that total ener- 
gies improve monotonically with decreasing e (or increas- 
ing Hilbert space size M(e)). Use can be made of 1/M(e) 
to extrapolate to the e — > limit. Results of such linear, 
two-point extrapolations for the N = 64 and 128 cases 
are included in Table || and generally improve the results 
from the largest value of M(e). Note, however, that the 
extrapolations are not variational in general, nor are any 
of the results for the gaps, which are the difference of two 
total energies. 



IV. USE OF THE TRANSLATION OPERATOR 

A feature of the four-block method is that the reduced 
Bloch symmetry of the four-block ring can be used to 
explicitly target states in momentum sectors other than 
k = and k = ir. For example, by constructing su- 
perblock states with a phase of —1 under shifts: 

\ni) <g> \n 2 ) (g> |n 3 ) <g> |n 4 )(g) — > \n 2 ) ® |n 3 ) ® |n 4 ) ® \ni), 

the k — 2ir/N symmetry sector can be targeted directly. 

For example, we consider the exactly solvable, 
dimerised XY spin chain: 

N 

H = 2j2(S?Sf +1 + SfSf +1 ) (1 + (-1)'**) , (14) 

i=l 

where 8 is the dimerisation parameter. The exact ex- 
citation spectrum in any jaaomcntum sector is readily 
obtained for finite latticesca. Table || shows the con- 
vergence of the four-block method for the energy gap 
Att/16, from the ground state to the lowest excitation 
in the k = it/16 sector, for the TV = 64 site ring with 
S = 0.2. Note that A^ng is resolved to within around 
0.01%. The reduced Bloch symmetry may prove to be 
useful in mapping out the excitation spectrum for, say, a 
64 site lattice. 
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V. APPLICATION OF THE FOUR-BLOCK 
METHOD TO SYSTEMS WITH A LARGE 
NUMBER OF DEGREES OF FREEDOM PER 
SITE 

In the above sections we have shown that the four- 
block method can provide accurate determinations of en- 
ergy gaps in spin models on large, periodic rings. How- 
ever, the real utility in the method lies in its ability to 
deal with systems with a large number of degrees of free- 
dom per site. To illustrate this we consider the XY spin 
chain interacting with dispersionless quantum phononstj. 
Here we are mainly concerned with demonstrating the 
convergence of the four-block method, rather than per- 
forming a comprehensive study of the model. The Hamil- 
tonian is given by 

n = J2{ 1 + y{ b Ui + b ^- b t- b *) 

i 

x (Sf^+^+^iK (15) 

i 

where bi destroys a phonon of frequency to and g is 
the spin-phonon coupling. This model has been stud- 
ied by ffron and Moukouri using the standard DMRG 
methods. The model undergoes a Kosterlitz-Thouless 
(K-T) transition at some critical coupling g c from a Lut- 
tinger liquid phase (g < g c ) with gapless excitations to a 
gapped, dimerised piiase (g > g c ) with a doubly degener- 
ate ground state. Ir£3 the gap A in the dimerised region 
is determined as a function of g in the thermodynamic 
limit by performing DMRG calculations cm large lattices. 
This data is fitted to Baxter's K-T formEd 



A = 



Vi 2 



9o 



: exp 



Vi 2 



fJc 



(16) 



in order to determine g c . 

Here we show that the four-block method can be used 
to accurately calculate a number of energy gaps on fi- 
nite, periodic rings. The examination of the crossover of 
these gaps allows an accurate determinatioia^i f . th e crit- 
ical point for models with a K-T transitionLaOCa. The 
good quantum numbers for this model that can be ex- 
ploited by the four-block method are the total z spin, 
(Sj^, and, in the = sector, the spin-flip symmetry 
(13|). In addition, the reduced Bloch symmetry and re- 
flection symmetry can be used. In addition to the ground 
state energy Eqs = E (S^ = 0, T = 1), three gaps are 
considered. In the notation of Nomura, these are: the 
doublet gap: A doub i c t = E (S% = ±1) - E GS ; the dimer 
gap: A dimor = J5 x (5f = 0,f = 1) - £ G a^and the Neel 



gap: 



iNcol 



E (S$ = 0,T = -1) 



E, 



GS>- 



The first stage of the four-block method requires the 
exact diagonalisation of a four-site ring. Convergence 
results for the various gaps in the four-site system are 
given in Table III for the ui = 10, g = 2.4 case. Note 



that in this case convergence with m, the number of 
bare phonon modes retained per site, is very rapid. If 
u> (the energy needed to create a phonon excitation) is 
decreased and/or the coupling g is increased, the value of 
m required for convergence will increase. In calculations 
performed on various electron-phonon modelscSEJ the re- 
quired value has not proved prohibitive i.e. m < 30. If 
this stage of the calculation does present a problem then 
global Hilbert space truncation can be used, i.e., rather 
than placing a limit m on the phonon number for each 
site, the sum of the phonon numbers from all sites is re- 
stricted. This simple step dramatically reduces the size 
of the Hilbert space required for convergence. Failing 
this, the local Hilbert space reduction method of White 
and co-workerstJ can be used to build an efficient basis 
for the f-site ring by starting with a two-site system and 
forming a single-site density matrix. Even the two-site 
calculation can be made more efficient by using a coher- 
ent state basis rather than simply bare phonon states. 

The convergence with e of subsequent stages of the 
four-block method — N = 8, 16 and 32 — is illustrated in 



Table [V. The convergence is sufficiently rapid that the 



data can be used in finite-size scaling studies (the gaps 
in the N = 32 case are resolved to within around 0.01%). 
As an example we consider the critical coupling g c . g c 
can be obtained as the limiting value of g c (N), where 
g c (N) denotes the finite-size crossover from the gapless 
spin-fluid phase to the dimerised phase, and is fixed by 
the conditionE£l 

^doublet — ^dimer- 

(17) 

That is, the lowest excitation in the fluid phase is the dou- 
blet, whose energy gap vanishes in the bulk limit, whilst 
in the dimer phase the dimer excitation becomes degen- 
erate with the ground state in the bulk limit, whereas 
the doublet energy gap. approaches a non-zero limit, the 
energy gap A used inE3. 

Plots of Adoubict — Adimor versus g in the u = 10 case 
are given in Fig. || for various values of N. A simple 
quadratic fit of the data in Fig. || gives estimates for 
g c (N) which are tabulated in Table. W. Note the very 
rapid convergence of g c (N) with TVljT^ For smaller val- 
ues of uj there is strongjmixing between fermion-like and 
phonon-like excitationsEj in the dimer and Neel sectors 
and it is only for large lattices that g c (N) converged^, 
when the characteristic electron gap, 2tt/N, falls below 
w, the energy required to create a phonon excitation. 

Taking into account the discretisation and fitting er- 
rors, we can safely estimate g c = 2.41(3). This result 
is to be compared with the.. result g c ~ 2.9, obtained 
from the phase boundary ir£j. The discrepancy between 
the two results is probably due to the problematic na- 
ture of fitting the infinite system gap A to (|l6|). That 
is, three parameters, a, b and g c , must be obtained from 
the non- linear fit, and it is very difficult to determine A 
accurately near g = g c . This is because A is extremely 
small for values of g even substantially higher than <? c , 
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due to the essential singularity in (|T^). Determining 
such small gaps from finite-size scaling is very difficult 
as very large lattices are required in order to observe the 
crossover from the initial algebraic scaling with N to the 
exponential scaling expected for gapped systems. This 
is further complicated by the fact that open^rather than 
periodic, boundary conditions were used inEa and by the 
presence of substantial DMRG truncation error for long 
chains and small gaps. These factors could well lead to 
an overestimation of g c . 

Finally, we perform a consistency check on our hypotb-. 



esis 
an> 



rat the transition is of the K-T type. Following 
we define parameters vq and v\ according to 
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tt 2 vq 
6N 



(18) 



and 



[2A d 



ablet 



iNool 



TTVi 
If 



(19) 



where is the bulk ground state energy per site, vq 
and V\ are extracted from the finite-size scaling of Eqs 
and the gaps. Two-point extrapolations of the TV = 16 
and N = 32 results give v /vi = 0.9980 for g = 2.4 sa g c . 
This is highly consistent with the result vq/vi p,l which 
should hold at the fluid-dimer transition point£3. 



VI. SUMMARY 



In this paper a new variant of White's density matrix 
renormalisation group (DMRG) method was presented. 
The algorithm was designed for the performance of finite- 
size scaling studies of one-dimensional quantum lattice 
models with a large number of degrees of freedom per site. 
We call the technique the "four-block method" , because 
superblocks consisting of four identical blocks are used as 
opposed to the standard DMRG method which uses two. 
The four-block method was shown to recover exact Bethe 
ansatz results for spin-1/2 Heisenberg rings of up to 128 
sites with good accuracy. It was shown that partial use 
can be made of the Bloch symmetry so that momentum 
sectors other than k = and k = it can be targeted 
directly. 

However, the real utility of the four-block method lies 
in its ability to treat systems with a large number of de- 
grees of freedom per site such as electron-phonon models. 
This was demonstrated by applying the method to the 
XY spin-Peierls model. It was shown that the method 
accurately resolves a number of finite-system energy gaps 
in this model, and, using finite-size scaling, the critical 
coupling was accurately determined for one particular 
phonon frequency. In future studies-,the Heisenberg spin- 
Peierls model will be investigated!^!. Extensions of the 
four-block method to higher dimensions are also being 
pursued. 
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TABLE I. The ground state energy, -Egs> an d the singlet 
and triplet gaps, A ss and A st , of N = 32, 64 and 128 site 
spin- 1/2 antiferromagnetic Heisenberg rings calculated using 
the four-block method for a number of cutoff parameters e. 
The size of the superblopk Hilbert space is M(e). Exact re- 
sults are taken from ref.EJ. Two-point extrapolations to the 
e — > limit (using 1/M(e)) are included for the N = 64 and 
N = 128 cases. 
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TABLE II. The energy gap A^/ 16 from the ground state 
to the first excited state in the k = -n/ 16 momentum sector 
for the N = 64 site dimerised XY model with dimerisation 
S — 0.2 calculated using the four-block method for a number 
of cutoff parameters e. The size of the superblock Hilbert 
space is M(e). 



e M(e) A^ /16 

10 _i2 2055402 0.90596 

10" 13 3882712 0.90476 

10" 14 6984756 0.90429 

10" 15 12074567 0.90411 

Exact 0.90401 



TABLE III. Convergence of various energy gaps at the 
first stage (exact diagonalisation of a four-site ring) of the 
four-block method for the XY spin-Peierls model with ui — 10 
and g = 2 AO. m is the number of bare phonon levels retained 
per site in the truncated basis and M(m) = 6m 4 is the size 
of the four-site Hilbert space. 
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TABLE IV. Convergence of various energy gaps calculated 
using the four-block method for the XY spin-Peierls model 
with to — 10 and g — 2.40 and TV = 8, 16 and 32 sites. Here e 
is the cutoff parameter which determines the accuracy of the 
method and M(e) is the size of the superblock Hilbert space. 
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TABLE V. Convergence with lattice size N of the crossover 



coupling g c (N) for the XY spin-Peierls model in the w = 10 
case. 
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9c(N) 
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2.5308 
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2.4338 


16 


2.4049 


32 


2.4083 
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(b) Stage 2 




(c) Stage 3 

FIG. 1. Schematic illustration of the four-block method, 
(a) The first step is to diagonalise a four-site ring. A reduced 
density matrix is defined for a two-site subsystem, (b) The 
density matrix eigenstates are then used as a truncated basis 
for the two-site blocks used in the second iteration, which 
consists of four t-site blocks, (c) Again, a reduced density 
matrix is defined for a two-block (four-site) subsystem and 
the density matrix eigenstates are used to form a truncated 
basis for the four-site blocks used in the third iteration. The 
size of the lattice (superblock) doubles at each iteration. 
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FIG. 2. The difference of the doublet and dimer gaps as a 
function of g for the XY spin-Peierls model with uj — 10 as 
determined by the four-block method for N = 4 (dotted line) , 
N = 8 (dot-dashed line), N = 16 (dashed line), and N = 32 
(full line). The crossover point g c (N) is determined by the 
condition Adoubiet = Adimcr (the intercepts of the curves with 
the horizontal axis) and converges rapidly to the critical point 
g c as N — > oo. 
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